/**
 * visTopoClustering.cpp
 * ubuntu8.04,g++4.3.2
 * 18/3/2009,VISUALIZATION_CGCAD_THSS_THU_BJ
 * Author : Xinlai,Lu
 *
 * Modification History:
 *
 */

/********@@@@@@@@@@@ PAY ATTENTION !!! @@@@@@@@@@@@**********/

#define GEN_STREAMLINES
//#define COPY_STREAMLINES
#define CLUSTERING_HI

/********@@@@@@@@@@@ PAY ATTENTION !!! @@@@@@@@@@@@**********/

#include "visTopoClustering.h"

namespace visualization
{
void visTopoClustering::parseParameter(string param)
{
	ifstream paramFile;
	paramFile.open(param.c_str());
	if(!paramFile)
	{
		cout << "Fail to open the specified parameter file, use the default one instead!" << endl;
		paramFile.open("default.param");
	}
	string parameters;
	while(getline(paramFile, parameters))
	{
		if('#' == parameters[0])continue;
		string data;
		size_t pos = parameters.find(' ');
		data = parameters.substr(pos + 1, parameters.length());
		parameters = parameters.substr(0, pos);
		double vlu = atof(data.c_str());
		if("m_preClusteringThreshold" == parameters) m_preClusteringThreshold = vlu;
		if("m_ClusteringThreshold" == parameters) m_ClusteringThreshold = vlu;
		if("m_CenterNum" == parameters)	m_CenterNum = size_t(vlu);
		if("m_Iterations" == parameters) m_Iterations = size_t(vlu);
		if("m_Step" == parameters) m_Step = vlu;
		if("m_MaxTime" == parameters) m_MaxTime = size_t(vlu);
		if("m_Segments" == parameters) m_Segments = size_t(vlu);
		if("m_Scale" == parameters) m_Scale = vlu;
		if("m_SourcesNum" == parameters) m_SourcesNum = size_t(vlu);
		if("m_Interval" == parameters) m_Interval = size_t(vlu);
		if("m_Threshold" == parameters) m_Threshold = vlu;
		if("m_UseThreshold" == parameters) m_UseThreshold = size_t(vlu);
		if("m_DuplicateThshold" == parameters) m_DuplicateThshold = vlu;
		if("m_UseRemoveDuplicate" == parameters) m_UseRemoveDuplicate = size_t(vlu);
		if("m_DistanceThshold" == parameters) m_DistanceThshold = vlu;
		if("m_Diavation" == parameters) m_Diavation = vlu;
		if("m_GridLen" == parameters) m_GridLen = vlu;
		if("m_SelfIncrementCenters" == parameters) m_SelfIncrementCenters = size_t(vlu);
		if("m_CrossIncrementCenters" == parameters) m_CrossIncrementCenters = size_t(vlu);
	}
	cout << setiosflags(std::ios::left);

	cout << setw(25) << "Parameters	:	" << endl;
	cout << setw(25) << "m_ClusteringThreshold	:	" << m_ClusteringThreshold << endl;
	cout << setw(25) << "m_preClusteringThreshold	:	" << m_preClusteringThreshold << endl;
	cout << setw(25) << "m_CenterNum : " << m_CenterNum << endl;
	cout << setw(25) << "m_Iterations : " << m_Iterations << endl;
	cout << setw(25) << "m_Step : " << m_Step << endl;
	cout << setw(25) << "m_MaxTime : " << m_MaxTime << endl;
	cout << setw(25) << "m_Segments : " << m_Segments << endl;
	cout << setw(25) << "m_Scale : " << m_Scale << endl;
	cout << setw(25) << "m_SourcesNum : " << m_SourcesNum << endl;
	cout << setw(25) << "m_Interval : " << m_Interval << endl;
	cout << setw(25) << "m_Threshold : " << m_Threshold << endl;
	cout << setw(25) << "m_UseThreshold : " << m_UseThreshold << endl;
	cout << setw(25) << "m_DuplicateThshold : " << m_DuplicateThshold << endl;
	cout << setw(25) << "m_UseRemoveDuplicate : " << m_UseRemoveDuplicate << endl;
	cout << setw(25) << "m_DistanceThshold : " << m_DistanceThshold << endl;
	cout << setw(25) << "m_Diavation : " << m_Diavation << endl;
	cout << setw(25) << "m_GridLen : " << m_GridLen << endl;
	cout << resetiosflags(std::ios::left);

	paramFile.close();
	cout << "Input the parameter file for streamline generation..." << endl;
	string pfile;
	cin >> pfile;
	m_ClusteringMethod.setParamFile(pfile);
}

void visTopoClustering::clusterDomain()
{
	cout << endl << "Begin clustering in domain..." << endl;
	cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl << endl;
	UCDInputFile file(m_FileName);
	m_Mdl = new Model;
	file.FGetModel(*m_Mdl);
	m_Mdl->buildKDTree();

#ifdef SCANLINE_FILL
	m_Mdl -> getBound(xmin, xmax, ymin, ymax, zmin, zmax);
	size_t dataPoint = m_Mdl -> getPntNumberOfModel();
	dataPoint = size_t(sqrt((double)dataPoint));
	xres = yres = size_t((dataPoint + 1) * 2.5);
	m_xintv = (xmax - xmin)/xres;
	m_yintv = (ymax - ymin)/yres;
	m_pointNum = xres * yres;	
#endif

	setDataField(m_Mdl);
	topoStructure();
	buildKDTree();
	roamDomain();
	preProcessSources();
	cout << "Sources Number now	:	" << endl;
	for(size_t i = 0; i < saddles.size(); i++)
	{
		cout << "Saddle	" << i << " : " << endl;
		for(size_t j = 0; j < QUADRANTS; j++)
		{
			cout  << m_SourcesInDomains[i].sources[j].size() << ",	";
		}
		cout << endl;
	}

	m_ClusteringMethod.setDataField(m_Mdl);
	m_ClusteringMethod.setParameters(m_Step, m_MaxTime);
	m_ClusteringMethod.setSegment(m_Segments);
	
	if(!m_ClusterCenters) 
	{
		cout << "Memory out at line " << __LINE__ << " of file " << __FILE__ << endl;
		cout << "Stop!" << endl;
		exit(-1);
	}

	time_t begin, end;
	begin = clock();

#ifdef GEN_STREAMLINES
	for(size_t i = 0; i < saddles.size(); i++)
	{
		cout << endl << "Saddle	" << i << " :" << endl;
		for(size_t j = 0; j < QUADRANTS; j++)
		{
			if(0 == m_SourcesInDomains[i].sources[j].size())
			{
				m_ClusterCenters[i * QUADRANTS + j].center = NULL;
				cout << "Empty sources... ... ";
				continue;
			}
			cout << j << "... ...  ";

			m_ClusteringMethod.setSources(m_SourcesInDomains[i].sources[j]);
			m_ClusteringMethod.computeStreamlinesInDomain();

#ifdef COPY_STREAMLINES
			vector<Point> tmp;
		  for(size_t m = 0; m < m_ClusteringMethod.getStreamlinesNumber(); m++)
		  {
		  	tmp.clear();
		  	for(size_t n = 0; n < m_ClusteringMethod.getStreamlineSize(m); n++)
		   	{
		   		tmp.push_back(m_ClusteringMethod.getPointOnStreamline(m, n));
		   	}
		  	m_Streamlines.push_back(tmp);
		  }
#endif
	
#ifdef CLUSTERING_HI
//			m_ClusteringMethod.computeLength();
			m_ClusteringMethod.initialize();
			m_ClusteringMethod.partition(m_ClusterCenters[i*QUADRANTS+j].cntNum,m_preClusteringThreshold,m_ClusteringThreshold);
			m_ClusterCenters[i * QUADRANTS + j].center = m_ClusteringMethod.getCenters();
#endif
		}
	}
#endif
	end = clock();
	cout << endl << "TIME USED FOR CLUSTERING	:	" << (double)(end - begin)/CLOCKS_PER_SEC << endl;
}

void visTopoClustering::topoStructure()
{
	cout << "Computing topology structure..." << endl;
	cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl << endl;
	setRdAndSmpRate(1, 4);
	setZero(1e-6);
	setDelt(1e-3);
	setDiavationFromSaddle(1);
	setNeighbourNum(10);
	setDisError(4);
	determineCriticalPoint();
	preProcessSaddles();
	integrateCurves(m_Forward, m_Backward);
	dataField->getBound(m_XMin, m_XMax, m_YMin, m_YMax, m_ZMin, m_ZMax);
}
void visTopoClustering::buildKDTree()
{
	cout << endl << "Building kd-tree..." << endl;
	cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl << endl;
	/* Since the data field is two dimensional, there're four boundaries. To store
	 * the boundaries, each of them is sampled by boundaryPntNum points.
	 */
	m_Saddles = saddles.size();
	const size_t boundaryPntNum = 2000;
	const size_t curveSegments = 2000;
	size_t curvePntNum = m_Saddles * 4 * (curveSegments + 1);
	/**The integration curves do not start exactly from the saddle points but with a 
	 * diavation. Therefor, the segments linking the source pairs are needed to build
	 * a kd-tree of the whole boundary.
	 */
	const size_t curveSourcesLink = 500;
	// The number of points passed to the kdTree must be exactly right!
	//intgCurvesAndBoundaries = new kdTree(curvePntNum + boundaryPntNum*4 + curveSourcesLink*m_Saddles*2);
	intgCurvesAndBoundaries = new kdTree(curvePntNum + boundaryPntNum*4 + curveSourcesLink*m_Saddles*2);

	// (m_XMin, m_YMin) ---> (m_XMax, m_YMin)
	// (m_XMin, m_YMax) ---> (m_XMax, m_YMax)
	// (m_XMin, m_YMin) ---> (m_XMin, m_YMax)
	// (m_XMax, m_YMin) ---> (m_XMax, m_YMax)
	for(size_t i = 0; i < boundaryPntNum; i++)
	{
		double x = (m_XMax - m_XMin)/(boundaryPntNum-1)*i;
		double y = (m_YMax - m_YMin)/(boundaryPntNum-1)*i;
		intgCurvesAndBoundaries->store(x, m_YMin, 0);
		intgCurvesAndBoundaries->store(x, m_YMax, 0);
		intgCurvesAndBoundaries->store(m_XMin, y, 0);
		intgCurvesAndBoundaries->store(m_XMax, y, 0);
	}

	visStreamlineClustering toResampleIntgCurves;
	for(size_t i = 0; i < m_Saddles * 2; i++)
	{
		vector<double> len1, len2;
		Point pnt1, pnt2;
		toResampleIntgCurves.computeLength((*m_Forward)[i], len1);
		toResampleIntgCurves.computeLength((*m_Backward)[i], len2);
		for(size_t seg = 0; seg <= curveSegments; seg++)
		{
			toResampleIntgCurves.getPointOnStreamline((*m_Forward)[i], len1, curveSegments, seg, pnt1);
			toResampleIntgCurves.getPointOnStreamline((*m_Backward)[i], len2, curveSegments, seg, pnt2);
			intgCurvesAndBoundaries->store(pnt1.getX(), pnt1.getY(), 0);
			intgCurvesAndBoundaries->store(pnt2.getX(), pnt2.getY(), 0);
		}

		if(0 == i % 2)
		{
			Point src1 = ForwardSources[i]; 
			Point src2 = ForwardSources[i+1] - src1;
			src2 /= curveSourcesLink;
			Point src3 = BackwardSources[i]; 
			Point src4 = BackwardSources[i+1] - src3;
			src4 /= curveSourcesLink;
			for(size_t k = 0; k < curveSourcesLink; k++)
			{
				src1 += src2;
				src3 += src4;
				intgCurvesAndBoundaries->store(src1.getX(), src1.getY(), src1.getZ());
				intgCurvesAndBoundaries->store(src3.getX(), src3.getY(), src3.getZ());
			}
		}
	}
	intgCurvesAndBoundaries->treeBalance();
}
void visTopoClustering::roamDomain()
{
	cout << endl << "Roaming the domains to determine streamline seeds..." << endl;
	cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl << endl;
	m_SourcesInDomains = new SourcesAroundSaddle[m_Saddles];
	m_ClusterCenters = new ClusterCenters[m_Saddles * QUADRANTS];
	/**
	 *source2 v1    source1
	 *   \    |   /
	 *     \  | /
	 * -------|------ v2
	 *      / | \
	 *    /   |   \
	 *source3    source4
	 */
	size_t index = 0;
	for(size_t sdl = 0; sdl < m_Saddles; sdl++)
	{
		Point source1(saddles[sdl].x + (saddles[sdl].eigenvec[0] + saddles[sdl].eigenvec[2])*m_Scale,
								  saddles[sdl].y + (saddles[sdl].eigenvec[1] + saddles[sdl].eigenvec[3])*m_Scale);
		Point source2(saddles[sdl].x + (saddles[sdl].eigenvec[0] - saddles[sdl].eigenvec[2])*m_Scale,
								  saddles[sdl].y + (saddles[sdl].eigenvec[1] - saddles[sdl].eigenvec[3])*m_Scale);
		Point source3(saddles[sdl].x + (- saddles[sdl].eigenvec[0] - saddles[sdl].eigenvec[2])*m_Scale,
								  saddles[sdl].y + (- saddles[sdl].eigenvec[1] - saddles[sdl].eigenvec[3])*m_Scale);
		Point source4(saddles[sdl].x + (- saddles[sdl].eigenvec[0] + saddles[sdl].eigenvec[2])*m_Scale,
								  saddles[sdl].y + (- saddles[sdl].eigenvec[1] + saddles[sdl].eigenvec[3])*m_Scale);

		m_GridInterval = 0;
		sourcesInDomain(m_SourcesInDomains[sdl].sources[0], source1, 'O');

		m_GridInterval = 0;
		sourcesInDomain(m_SourcesInDomains[sdl].sources[1], source2, 'O');

		m_GridInterval = 0;
		sourcesInDomain(m_SourcesInDomains[sdl].sources[2], source3, 'O');

		m_GridInterval = 0;
		sourcesInDomain(m_SourcesInDomains[sdl].sources[3], source4, 'O');
		for(size_t m = 0; m < QUADRANTS; m++)
		{
			m_ClusterCenters[sdl * QUADRANTS + m].cntNum = m_CenterNum;
		}
	}
}

void visTopoClustering::preProcessSources()
{
	size_t toDelSdl, toDelQd;
	for(size_t i = 0; i < m_Saddles; i++)
	{
		for(size_t j = 0; j < QUADRANTS; j++)
		{
			for(size_t s = 0; s < m_SourcesInDomains[i].sources[j].size(); s++)
			{
				bool del = false;
				// 
				for(size_t ii = 0; ii < m_Saddles; ii++)
				{
					if( del ) break;
					for(size_t jj = 0; jj < QUADRANTS; jj++)
					{
						if( del ) break;
						if( i == ii && j == jj ) continue;
						for(size_t ss = 0; ss < m_SourcesInDomains[ii].sources[jj].size(); ss++)
						{
							Point pnt1 = m_SourcesInDomains[i].sources[j][s];
							Point pnt2 = m_SourcesInDomains[ii].sources[jj][ss];
							double dist = pnt1.Distance( pnt2 );
							if( dist < m_GridLen * 0.9 )
							{
								del = true;
								toDelSdl = m_SourcesInDomains[i].sources[j].size() > m_SourcesInDomains[ii].sources[jj].size() ?
													 ii : i;
								toDelQd = m_SourcesInDomains[i].sources[j].size() > m_SourcesInDomains[ii].sources[jj].size() ?
													jj : j;
								if(i == ii)
								{
									m_ClusterCenters[ii*QUADRANTS+jj].cntNum += m_ClusterCenters[ii*QUADRANTS+jj].cntNum +
																														  m_SelfIncrementCenters;
								}
								else
								{	m_ClusterCenters[ii*QUADRANTS+jj].cntNum += m_CrossIncrementCenters; }
								break;
							}
						}
					}
				}
				if( del )	
				{
					Point tmp = m_SourcesInDomains[toDelSdl].sources[toDelQd][0];
					m_SourcesInDomains[toDelSdl].sources[toDelQd].clear();
//					m_SourcesInDomains[toDelSdl].sources[toDelQd].push_back(tmp);
//					m_ClusterCenters[toDelSdl*QUADRANTS+toDelQd].cntNum = 1;
				}
			}
		}
	}
	for(size_t i = 0; i < m_Saddles ; i++)
	{
		for(size_t j = 0; j < QUADRANTS; j++)
		{
			size_t sz = m_SourcesInDomains[i].sources[j].size();
			if(sz < m_ClusterCenters[i*QUADRANTS+j].cntNum)
			{
				m_ClusterCenters[i*QUADRANTS+j].cntNum = sz;
			}
		}
	}
}
//
//void visTopoClustering::postProcess()
//{
//	cout << endl << "Post-processing..." << endl;
//	const size_t nearestNum = 1;
//	const double disError = pow(m_XMin-m_XMax,2)+pow(m_YMin-m_YMax,2)+pow(m_ZMin-m_ZMax,2);
//	
//	visStreamlineClustering *computeLen = new visStreamlineClustering;
//	vector<double> len;
//	// Remove the cluster-center if it's very near to the boundary.
//	for(size_t i = 0; i < m_Saddles * QUADRANTS; i++)
//	{
//		if(!m_ClusterCenters[i].center) continue;
//		for(size_t k = 0; k < m_ClusterCenters[i].cntNum; k++)
//		{
//			double dis = 0;
//			computeLen -> computeLength( m_ClusterCenters[i].center[k + m_Iterations * m_ClusterCenters[i].cntNum], len );
//			for(size_t j = 0; j <= m_Segments; j++)
//			{
//				Point tmp;
//				computeLen->getPointOnStreamline( m_ClusterCenters[i].center[k + m_Iterations * m_ClusterCenters[i].cntNum], 
//																		      len, m_Segments, j, tmp);
//				nNearestNodes nNN(nearestNum);
//				nNN.setDistError(disError);
//				nNN.setSearchPnt(tmp.getX(), tmp.getY(), tmp.getZ());
//				intgCurvesAndBoundaries->locateNodes(&nNN, 1);
//				if(0 == nNN.found)
//				{
//					cout << "Warning !!!!!! " << endl;
//					continue;
//				}
//				dis += sqrt(nNN.dist2[1]);
//			}	
//			if(0 == len.size())
//			{
//				cout << "Warning!!! The length is zero!!!" << endl;
//				continue;
//			}
//			dis /= (m_Segments + 1);
//			if(m_Threshold > dis)
//			{
//				m_ClusterCenters[i].center[k + m_Iterations * m_ClusterCenters[i].cntNum].clear();
//			}
//		}
//	}
//	delete computeLen;
//}
//
//void visTopoClustering::removeDuplicate()
//{
//	cout << "Removing streamlines that are similar..." << endl;
//	visStreamlineClustering *computeSim = new visStreamlineClustering;
//	computeSim->setSegment(m_Segments);
//
//	const size_t times = m_CenterNum - 1;
//	size_t t = 0;
//	while(t++ < times)
//	{	
//		for(size_t i = 0; i < m_Saddles * QUADRANTS; i++)
//		{
//			if(!m_ClusterCenters[i].center) continue;
//			for(size_t m = 0; m < m_ClusterCenters[i].cntNum; m++)
//			{
//				size_t toRmQ, toRmCnt;
//				double sim = DBL_MAX;
//				for(size_t j = i; j < m_Saddles * QUADRANTS; j++)
//				{
//					if(!m_ClusterCenters[j].center) continue;		
//					for(size_t n = 0; n < m_ClusterCenters[j].cntNum; n++)
//					{
//						double tmp = computeSim->EuDistance( m_ClusterCenters[i].center[m+m_Iterations*m_ClusterCenters[i].cntNum],
//																								 m_ClusterCenters[j].center[n+m_Iterations*m_ClusterCenters[j].cntNum]
//																							 );
//						if(sim > tmp && !(i == j && m == n))
//						{
//							toRmQ = j;
//							toRmCnt = n;
//							sim = tmp;
//						}
//					}
//				}
//				if(m_DuplicateThshold > sim/*&& n != m_CenterNum - 1*/)
//				{
//					m_ClusterCenters[toRmQ].center[toRmCnt + m_Iterations * m_ClusterCenters[toRmQ].cntNum].clear();
//				}
//			}
//		}
//	}
//	delete computeSim;
//}

SOURCES visTopoClustering::sourcesInDomain(Point &sourceSeed)
{
	Point currentPos = sourceSeed;
	SOURCES srcs;
	srcs.push_back(currentPos);

	static size_t total = 0;
	size_t record = 0;

	const size_t nearestNum = 1;
	// To garant that a point will be found every time, the disError should be the length of the
	// diagnal of the data field.
	const double disError = pow(m_XMin-m_XMax,2)+pow(m_YMin-m_YMax,2)+pow(m_ZMin-m_ZMax,2);
	Point cnt;
	while(record < m_SourcesNum)
	{
		/* 1. Search the shortest distance "shortest" from the current position to the point "p" 
		 * 		in the kdtree.
		 * 2. Shitf the current postion by a distance "shortest" to the opposite direction of
		 * 		the vector pointing from "p" to the current position, and save the new position.
		 */
		record++;
		nNearestNodes nNN(nearestNum);
		nNN.setDistError(disError);
		nNN.setSearchPnt(currentPos.getX(), currentPos.getY(), currentPos.getZ());
		intgCurvesAndBoundaries->locateNodes(&nNN, 1);
		if(0 == nNN.found)
		{
			total++;
			cout << "Warning! Fails to find the shortest distance." << endl;
			cout << "Do it brutely..." << endl;
			Point tmp(1, 1, 0);
			currentPos += tmp;
			cout << "Done!" << endl;
			continue;
//			break;
		}
		Point pntOnTree(nNN.index[nearestNum]->pos[0],nNN.index[nearestNum]->pos[1],nNN.index[nearestNum]->pos[2]);
//1		currentPos = currentPos * 2 - pntOnTree;
//
//		Point theVector = currentPos - pntOnTree;
//		const double pi = 3.14;
//		double rd = rand()/double(RAND_MAX) * pi * 2;
//		theVector.rotate(rd);
//		currentPos += theVector;
		
//3   
		cnt += currentPos;
		Point tmp = currentPos - (cnt/record);
		double m = sqrt( pow(tmp.getX(),2)+pow(tmp.getY(),2)+pow(tmp.getZ(),2));
		if(0 == m) { currentPos = currentPos * 2 - pntOnTree; }
		else
		{
			double dist = currentPos.Distance(pntOnTree);
			if(dist < 10) { currentPos = currentPos * 2 - pntOnTree;	}
			else
			{
				tmp *= dist/m;
				currentPos += tmp;	
			}
		}

	  if(0 == record%m_Interval)
		{	srcs.push_back(currentPos); }
	}

	return srcs;
}

void visTopoClustering::sourcesInDomain(SOURCES &srcs, Point sourceSeed, char direction)
{
	vector<Point>::iterator iter;
	iter = find(srcs.begin(), srcs.end(), sourceSeed);
	if(iter != srcs.end()) return;
	if(srcs.empty()) srcs.push_back(sourceSeed);

	const size_t nearestNum = 1;
	const double disError = ( pow(m_XMin-m_XMax,2)+pow(m_YMin-m_YMax,2)+pow(m_ZMin-m_ZMax,2) );
	nNearestNodes nNN(nearestNum);
	nNN.setDistError(disError);
	nNN.setSearchPnt(sourceSeed.getX(), sourceSeed.getY(), sourceSeed.getZ());
	intgCurvesAndBoundaries->locateNodes(&nNN, 1);
	if(0 == nNN.found)
	{
		cout << "Warning! Fails to find the shortest distance." << endl;
		cout << "Point	:	" << sourceSeed.getX() << "	" << sourceSeed.getY() << endl;
		return;	
	}

	double dist = sqrt(nNN.dist2[1]);
	if(dist < m_DistanceThshold){	return;	}
		
	srcs.push_back(sourceSeed);
	Point tmpX(m_GridLen, 0, 0), tmpY(0, m_GridLen, 0);
	sourcesInDomain(srcs, sourceSeed + tmpX, 'U');
	sourcesInDomain(srcs, sourceSeed - tmpX, 'U');
	sourcesInDomain(srcs, sourceSeed + tmpY, 'U');
	sourcesInDomain(srcs, sourceSeed - tmpY, 'U');
}

// new...
void visTopoClustering::postProcess()
{
	cout << endl << "Post-processing..." << endl;
	const size_t nearestNum = 1;
	const double disError = pow(m_XMin-m_XMax,2)+pow(m_YMin-m_YMax,2)+pow(m_ZMin-m_ZMax,2);
	
	visStreamlineClustering *computeLen = new visStreamlineClustering;
	vector<double> len;
	// Remove the cluster-center if it's very near to the boundary.
	for(size_t i = 0; i < m_Saddles * QUADRANTS; i++)
	{
		if(!m_ClusterCenters[i].center) continue;
		for(size_t k = 0; k < m_ClusterCenters[i].cntNum; k++)
		{
			double dis = 0;
			computeLen -> computeLength( m_ClusterCenters[i].center[k], len );
			for(size_t j = 0; j <= m_Segments; j++)
			{
				Point tmp;
				computeLen->getPointOnStreamline( m_ClusterCenters[i].center[k], 
																		      len, m_Segments, j, tmp);
				nNearestNodes nNN(nearestNum);
				nNN.setDistError(disError);
				nNN.setSearchPnt(tmp.getX(), tmp.getY(), tmp.getZ());
				intgCurvesAndBoundaries->locateNodes(&nNN, 1);
				if(0 == nNN.found)
				{
					cout << "Warning !!!!!! " << endl;
					continue;
				}
				dis += sqrt(nNN.dist2[1]);
			}	
			if(0 == len.size())
			{
				cout << "Warning!!! The length is zero!!!" << endl;
				continue;
			}
			dis /= (m_Segments + 1);
			if(m_Threshold > dis)
			{
				m_ClusterCenters[i].center[k].clear();
			}
		}
	}
	delete computeLen;
}

void visTopoClustering::removeDuplicate()
{
	cout << "Removing streamlines that are similar..." << endl;
	visStreamlineClustering *computeSim = new visStreamlineClustering;
	computeSim->setSegment(m_Segments);

	const size_t times = 1;
	size_t t = 0;
	while(t++ < times)
	{	
		for(size_t i = 0; i < m_Saddles * QUADRANTS; i++)
		{
			if(!m_ClusterCenters[i].center) continue;
			for(size_t m = 0; m < m_ClusterCenters[i].cntNum; m++)
			{
				size_t toRmQ, toRmCnt;
				double sim = DBL_MAX;
				for(size_t j = 0; j < m_Saddles * QUADRANTS; j++)
				{
					if(!m_ClusterCenters[j].center) continue;		
					for(size_t n = 0; n < m_ClusterCenters[j].cntNum; n++)
					{
						double tmp = computeSim->EuDistance( m_ClusterCenters[i].center[m],
																								 m_ClusterCenters[j].center[n]
																							 );
						if(sim > tmp && !(i == j && m == n))
						{
							if( 1 == m_ClusterCenters[j].cntNum )
							{
								toRmQ = i;
								toRmCnt = m;
							}
							else
							{
								toRmQ = j;
								toRmCnt = n;
							}
							sim = tmp;
						}
					}
				}
				if(m_DuplicateThshold > sim && toRmCnt != 0)
				{
					m_ClusterCenters[toRmQ].center[toRmCnt].clear();
				}
			}
		}
	}
	delete computeSim;
}

void visTopoClustering::makeArrow()
{
	cout << "Making arrows..." << endl;
	
	size_t sz = 0;
	for(size_t i = 0; i < m_Saddles * QUADRANTS; i++)
	{
		if(!m_ClusterCenters[i].center) continue;
		sz += m_ClusterCenters[i].cntNum;
	}
	m_Arrows = new Point[sz * 3];
	size_t count = 0;

	const double disError = ( pow(m_XMin-m_XMax,2)+pow(m_YMin-m_YMax,2)+pow(m_ZMin-m_ZMax,2) );
	visStreamlineClustering *computeCenter = new visStreamlineClustering;
	vector<double> *val = new vector<double>;
	for(size_t i = 0; i < m_Saddles * QUADRANTS; i++)
	{
		if(!m_ClusterCenters[i].center) continue;
		for(size_t j = 0; j < m_ClusterCenters[i].cntNum; j++)
		{			
			vector<double> len;
			computeCenter -> computeLength( m_ClusterCenters[i].center[j], len );	
			double length = len[ len.size() - 1 ] / 2;
			size_t index = computeCenter -> binarySearch( len, length );
			m_Arrows[count] = m_ClusterCenters[i].center[j][index];
			double pos[3] = {
				m_Arrows[count].getX(), 
				m_Arrows[count].getY(),
				m_Arrows[count].getZ()
			};
			m_Mdl->interpolatedValue(pos, 3, val, 4, disError);
			Point dir( (*val)[0], (*val)[1], (*val)[2] );
			m_Arrows[count].drawArrow( m_Arrows[count], dir, m_Arrows + count + 1);
			count += 3;
		}
	}
	delete computeCenter;
	delete val;
}

}//visualization

